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Abstract 

A detailed numerical study of two-dimensional flow past a circular cylinder at 
moderately low Reynolds numbers has been conducted using three different numerical 
algorithms for solving the time-dependent compressible Navier-Stokes equations. It 
was found that if the algorithm and associated boundary conditions were consistent 
and stable, then the major features of the unsteady wake were well-predicted. However, 
it was also found that even stable and consistent boundary conditions could introduce 
additional periodic phenomena reminiscent of the type seen in previous wind-tunnel 
experiments. However, these additional frequencies were eliminated by formulating 
the boundary conditions in terms of the characteristic variables. An analysis based on 
a simplified model provides an explanation for this behavior. 


‘Research was supported by the National Aeronautics and Space Administration under NASA Contract 
Nos NASl-18107 and NAS1-18605 while the first three authors were in residence at the Institute for 
Computer Applications in Science and Engineering (ICASE), NASA Langley Research Center, Hampton, 

VA 23665. 
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1. Introduction 

Historically, the unsteady wake generated by a circular cylinder in low-speed flow has 
been of great interest to computational fluid dynamicists as well as to theoretical and 
experimental aerodynamicists. The Reynolds-number range between 40 and 1000 has 
been of particular interest because it spans the transition from steady flow to unsteady 
wake flow dominated by the periodic shedding of vortices from the cylinder. The shedding 
frequency of these vortices increases with Reynolds number over this range, asymptotically 
approaching a constant value. However, Sreenivasan (1985) measured more than one 
distinct frequency in the shedding regime at low Reynolds numbers. In addition to the 
vortex-shedding frequency, he found clearly discernible, lower frequencies in the frequency 
spectrum for the streamwise velocity measured in the wake. These additional frequencies 
were not subharmonics of the primary shedding frequency. He concluded that this was a 
feature of the initial stages of transition to turbulence in agreement with the route to chaos 
described by Ruelle and Takens (1971). Sirovich (1985) suggested that these additional 
modes of oscillation could be described theoretically in terms of the classical von Karman 
vortex street. On the other hand, based on measurements of vibrating and nonvibratmg 
cylinders, van Atta and Gharib (1987) concluded that the additional frequencies found by 
Sreenivasan were due to the aeroelastic coupling of the vortex wake with cylinder vibration 
modes. For a nonvibrating cylinder, they found no spectral peaks other than the primary 
Strouhal shedding frequency. Recently, however, Sreenivasan (1989) determined that the 
cylinders did not vibrate in his previous experiments. 

Conflicting results have also come from the numerical studies of this flow. Karniadakis 
and Triantafyllou (1989) found no secondary frequency in their computation of the incom- 
pressible flow past a circular cylinder using the spectral-element method. Moreover, they 
were able to excite a secondary mode by introducing an external forcing function into the 
momentum equations. Townsend et al. (1987) did find a secondary frequency in their 
finite-difference computation of low-speed compressible flow past a circular cylinder. This 
low frequency was found in the frequency spectrum of the pressure at various points in 
the wake. However, the frequency was dependent upon the size of the solution domain, 
suggesting that it might have been a numerical effect. 

The present paper describes the results of further numerical studies which demonstrate 
that the secondary frequency found in the compressible-flow calculations was produced by 
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the far-field boundary conditions. Three different numerical methods were used to solve 
the time-dependent compressible Navier-Stokes equations. The first of these was the finite- 
difference technique used by Townsend et al. (1987), the second method was a fully-spectral 
method (Don 1989), while the third was a mixed spectral-finite-difference method (Don 
1989) which was a combination of the other two methods. Computations were made with 
the finite-difference code with two sets of far-field boundary conditions. These calculations, 
which were made for flow conditions in both the vortex-shedding regime and the steady- 
wake regime, demonstrate the effect of the boundary conditions on the frequency spectrum 
of the pressure in the wake. With the proper choice of far-field boundary conditions, 
no frequencies other than the primary shedding frequency and multiples of the primary 
frequency were found at a Reynolds number of 80. In addition, when the “improper” 
boundary conditions were used to compute the flow at a Reynolds number of 20, where the 
flow is known to be steady, then indeed no shedding frequency was predicted although the 
secondary frequency did appear in the solution. These findings were further substantiated 
by calculations with the highly-accurate spectral method. 

The details of the solution procedures, including boundary conditions, are described in 
§2 of the paper. All of the computed results are presented in §3. In §4, a mathematical 
analysis of boundary conditions for a model system is presented which shows how additional 
frequencies can be produced by the boundary conditions. Finally, a summary of the study 
and some concluding remarks are given in §5. 

2. Problem formulation and numerical methods 

2.1 Governing equations 

For the present compressible-flow calculations, the governing equations are those which 
describe the conservation of mass, momentum, and energy of an ideal fluid in the absence 
of external forces. The nondimensional form of these equations in a general curvilinear 
coordinate system is 

dU dF dG _ 1 (dF v dG v \ 

dt + d£ + dr] Re ref \ d£ + dr) ) ^ 

where 
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and 



dv dr^dv\ _ f d^du cty du\ 
\dyd^ + dyd rj ) \dxd£ dxdrj) ' 


The nondimensional equation of state is 


P = (7 - 1 )p t > 


( 2 . 6 ) 


(2.7) 


and the viscosity was given by the Sutherland law 

CiT § 

M “ c 2 + r 


( 2 . 8 ) 


where C x and C 2 are constants. The ratio of specific heats 7 was 1.4, and the Prandtl num- 
ber Pr was 0.72. These equations were nondimensionalized using the following reference 
quantities: U ref = U p re/ = Poo, Pr,/ = Pco^, and T ref = f/^/c u where 0*, and p*, are 
the free-stream velocity and density, respectively, and c v is the specific heat at constant 
volume. Therefore, the reference Reynolds number is defined as Re ref = U re fD/u ref where 
D is the diameter of the cylinder and u re f is the kinematic viscosity based on the reference 
temperature. However, it should be noted that the Reynolds number used in subsequent 
sections of the paper to characterize the flow conditions for which calculations were made 
is based upon free-stream quantities, i.e., Re oo X) — t^oo^/^oo- 


2.2 Grid 

A polar grid was used in the calculations. The grid was generated in the physical (x,y) 
plane and mapped numerically onto the computational (£,r)) plane. The £ coordinate 
corresponded to the circumferential direction and the rj coordinate to the radial direction. 
The grid was stretched in both directions so that grid points were clustered toward the 
cylinder in the radial direction and toward the wake region in the circumferential direction. 
For the finite-difference calculations, the basic grid had 122 points in the circumferential 
direction and 151 points normal to the body. The outer boundary of the grid was located 25 
diameters from the cylinder surface. For the computations with the spectral method, the 
grid had 64 and 48 points in the angular and radial directions, respectively. Calculations 
were made with the outer boundary at both 20 and 22.5 diameters away from the cylinder. 
For computations with the mixed spectral-finite-difference method, the grid had 70 and 
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J is the Jacobian of the transformation from Cartesian coordinates (x,y) to the general 
curvilinear coordinates (£, tj), 


and the specific internal 


J = 


dr) df dr\ 


dx dy dy dx 
energy, E, is defined by 


E = p 


T + 




( 2 . 2 ) 


( 2 . 3 ) 


The elements of the stress tensor are 



fd£ du dr) du\ fd£dv dr) dv\ 

dx dr] J \dyd£^ dy dr) J ’ 


( 2 . 4 ) 


fd£ du drj du\ fd£ dv dr) dv\ 
\dy d£ ^ dy dr) ) ^ \dx d£ ^ dx dt] ) ’ 


( 2 . 5 ) 
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100 points in the circumferential and radial directions, respectively. For this case, the far- 
field boundary was located 22.5 diameters from the cylinder surface. In the experiments 
of Sreenivasan (1985), the wind-tunnel walls were located 30 diameters above and below 
the cylinder 1 . 


2.3 Numerical methods 

The finite-difference method used in the present study was the original unsplit tech- 
nique of MacCormack (1969), which has second-order accuracy in both space and time. 
The method is an explicit, conditionally-stable, predictor-corrector scheme. Forward dif- 
ferences were used to approximate the derivatives of the fluxes in equation (2.1) in the 
predictor step, and backward differences were used in the corrector step. The first deriva- 
tives appearing in the viscous flux terms, F v and G u , were approximated with backward 
differences if the flux derivative was being approximated with a forward difference, and 
vice-versa. 

The spectral algorithm used both the Fourier and Chebyshev collocation methods. 
Because of the polar grid, the flow could be treated as being periodic at the boundaries 
of the solution domain in the circumferential direction. Thus, the Fourier collocation 
(pseudospectral) method was the natural choice for the circumferential direction. Since 
the flow in the radial direction was not periodic, the Chebyshev collocation method was 
used in that direction. Both of these methods were employed in the form of the matrix 
vector multiplication method instead of the more commonly used Fast Fourier Transform 
method. To improve the stability of the algorithm, a fourth-order exponential filter (Don 
1989) was applied to the differentiation and solution matrices. The solution was advanced 
in time using a second-order Runge-Kutta method. 

A third approach which combined these two methods was also used. In this case, the 
Fourier method was used in the circumferential direction and the finite-difference method 
in the radial direction. 

2.4 Boundary conditions 

A sketch of the solution domain in the physical plane is shown in figure 1. Periodic 
boundary conditions for the dependent variables were applied along the cut line in the 
1 Private communication. 
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wake for both the finite-difference and spectral methods. Standard solid-wall boundary 
conditions for viscous flow were applied at the cylinder surface. A no-slip condition was 
applied so that the velocity components, u and v, were specified to be zero at the surface. 
The wall temperature was held fixed at the value of the free-stream temperature. In 
the finite-difference code, the density at the wall was obtained from an extrapolation of 
interior-point values of pressure in the direction normal to the wall. In the spectral code, 
the density at the wall was computed as a part of the solution. 

The outer boundary is an arbitrarily-chosen boundary, introduced only to restrict the 
computational domain to a finite size. Since this is a boundary across which fluid passes 
either into or out of the computational region and since disturbances can propagate up- 
stream as well as downstream in a subsonic flow, careful consideration must be given to 
the specification of boundary conditions on this outer boundary. If this outer boundary is 
located sufficiently far from the cylinder surface, viscous effects are negligible in the flow 
crossing the boundaries except in the narrow region where the cylinder wake is located. 
As a result, the proper choice of boundary conditions can be found from an analysis of the 
inviscid form of the governing equations. 

For a subsonic inflow boundary, this analysis shows that there are three characteristics 
coming into the solution domain and one outgoing characteristic. Thus, there must be 
three quantities specified at an inflow boundary. Since there are four governing equations, 
the numerical method requires a fourth boundary condition in addition to the three re- 
quired for proper specification of the boundary conditions in a mathematical sense. In the 
finite-difference code, calculations were made with two different sets of inflow boundary 
conditions. The first set, which was used in the previous solutions reported by Townsend 
et al. (1987), was formulated in terms of the primitive variables u,v,p, and T. At the 
inflow boundary, free-stream values of the two velocity components and the temperature 
were specified, i.e., 


^inflow 

= Uoo 

(2.9) 

Vi n flow 

= v oo 

(2.10) 

Tin flow 

— T 

- 1 oo 

(2.11) 


The density was obtained from a zeroth-order extrapolation of pressure from the interior 
of the solution domain and the use of equation (2.7). In preliminary calculations in whic K 
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only the primary shedding frequency was of interest, this boundary condition gave results 
which were not significantly different from those obtained using an extrapolation of the 
outgoing characteristic to obtain pressure. 

At the outflow boundary, only one analytic boundary condition is needed since there is 
only one incoming characteristic. A typical choice is the specification of the static pressure. 
However, as demonstrated numerically by Rudy and Strikwerda (1980), this can cause 
waves to be reflected from the outflow boundary back into the solution domain, adversely 
affecting the solution in the interior of the domain. A nonreflecting boundary condition 
based on the work of Engquist and Majda (1977) was used in the present study. Thus, the 
pressure at the outflow boundary was found from a finite-difference approximation to the 
equation 

dp ( du dv\ 

a? ““a?) (212) 

where c is the nondimensional speed of sound given by 

c = y/iii-tyT . (2.13) 

This boundary condition was applied along the outflow boundary where the viscous wake 
crossed the boundary as well as the region immediately above and below the wake. Along 
the remainder of the outflow boundary, where the flow was essentially inviscid, the pressure 
was specified to be the free-stream pressure. The variables u,v, and T along the entire 
outflow boundary were obtained from zeroth-order extrapolation, and the density was then 
obtained from equation (2.7) using the boundary value of p and the extrapolated value of 
T. 

The second set of boundary conditions, which was also used in the spectral calculations, 
was based entirely on the characteristic variables. An analysis of the inviscid form of 
equation (2.1), linearized around the free-stream conditions, shows (Gottlieb, Lustman, 
and Streett 1984) that the characteristic variables, with corresponding eigenvalues a x = 


— *■ — # — * — * — *■ —* 
a 2 = U • N , a$ — U • N — c, and o 4 = U • N + c, are 

Ri = p^-(^pUoo-Uoo-M-Uoo + E^ (2.14) 

R 2 = (m u - PUoo)v v ~ - pv oo)*?* (2.15) 

R 3 = -{M-pU OQ )-N+ i 1 ^) (^pUoo-Uoo-M-Uoo + E^j (2.16) 

R, = (M — pUoo) • N + (^P^ 00 ' Uoo ~ M • Uoo + E^j (2.17) 
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where = (u^Uo,) is the free-stream velocity, M = (m u ,m„) = (pu,pv) is the local mo- 
mentum, Cqc is the free-stream speed of sound, and N — ( r} x , q y ) = (|^, |^)/ ^/(§f) 2 + (f^) 2 
is the unit outward normal vector. At the subsonic inflow boundary, the characteristic vari- 
ables corresponding to the three incoming characteristics were specified using free-stream 
values. These boundary values then became 

Ft = RiiPoo'&^E"), * = 1,2,3 (2.18) 

Furthermore, the required numerical boundary condition is 

Fl| Rii^Pnumt ^n»ni) (2.19) 

where values in F 4 are obtained by extrapolation from the interior of the domain in the 
finite-difference code and are obtained as a part of the solution in the spectral method. 
Similarly, at the outflow boundary, the characteristic variable corresponding to the incom- 
ing characteristic was specified, i.e., 

Fs — Rz (Poo 5 A^oo j Eqq ) (2.20) 

and the values for the other three characteristic variables were found numerically as at the 


inflow boundary, so that 

Ft — Ri{Pnumt Af num , E num ), * 1»2,4. (2-21) 

At each boundary, the conserved variables were then computed from the appropriate sys- 
tem of equations, giving 

p = lZl{F 1 + 4> 1 ) (2.22) 

L oo 

m u = (02 *7x "I" 03 Vy) “f” P^OO (2.23) 

m v = (02 'Hy - 03^) + pv 00 (2.24) 

1 w 

E = 0! - -pUw • Uoo + m u u oo + m v v <*> (2.25) 

z 

on the outer boundary, where 


<h = <h = \(F t -F s ), <h = F,. 
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3. Results and discussion 

Computations were made with all three methods for Mach 0.4 flow at a Reynolds 
number of 80. The initial flow field was specified to be free-stream flow everywhere except 
at the cylinder surface. Figure 2(a) shows pressure contours (isobars) from the finite- 
difference solution using the primitive-variable far-field boundary conditions at a point in 
time after the periodic vortex shedding had been fully established. The pressure field for 
finite-difference solution using the characteristic far-field boundary conditions is shown in 
figure 2(b). The corresponding solution for the fully-spectral code is given in figure 3. 

As the computations were made, the values of the pressure were saved at ten selected 
grid points in the flow field. At least four of these grid points were located in the wake 
region and the others were upstream of the cylinder. For the finite-difference code, the 
computed pressure data were saved every twenty tune steps, and for the fully-spectral 
and mixed spectral-finite-difference codes, because of the larger time step used in these 
methods, the data were saved after every time step. 

The set of pressure data at each point was analyzed for its frequency spectrum using 
a discrete Fast Fourier Transform (FFT) technique. First, the time-averaged pressure was 
subtracted from the pressure data to form a time sequence p, of N points, where N = 2 m 
for some positive integer m. Second, the Fourier coefficients F k of the transformation 

F k = for k = 0,1, . . . , N - 1 (3.1) 

j=0 

were obtained using a pre-coded FFT subroutine. A — 1 in the subroutine used to analyze 
the finite-difference results and A — 1 /N in the subroutine used to analyze the results 
from the fully-spectral and mixed-method codes. The power spectral density d k for wave 
number k was found then simply as 

d k = B\F k \ 2 where k = 0, 1 , ... ,N — 1. (3-2) 

B = 2/N for the finite- difference results and B = 1 for the others. The Strouhal number 
(nondimens ional frequency) corresponding to the vortex-shedding frequency was 

Sti = ki/iNAt), (3.3) 

where A t was the time step and k\ was the value of k such that 

d kl = max d k ■ 
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Similarly, St 2 is the Strouhal number corresponding to the secondary frequency. 

Note that, because the power density was found only for discrete frequencies, there 
was an inherent discretization error in determining the dominant frequencies in the flow. 
The actual peak frequency would be within half of the discretization, or ±l/(2iVAt}, of 
that given by the FFT. In some of the present cases, the resulting error in the Strouhal 
number could be as large as 0.003. In an effort to get a somewhat better estimate, for each 
dominant wave number a peak in the frequency spectrum was estimated using a three- 
point parabolic interpolation. These values were used in table 1, which summarizes the 
results found in the present study. 


Method 

R^oo t D 

Moo 

Domain 

Inflow BC 

Outflow BC 

Sh 

St% 

FD 

80 

0.4 

51D 

Primitive 

Primitive 

.1515 

.0165 

FD 

80 

0.2 

51D 

Primitive 

Primitive 

.1504 

.0347 

FD 

80 

0.1 

51D 

Primitive 

Primitive 

.1564 

.0721 

FD 

80 

0.4 

51D 

Characteristic 

Characteristic 

.1518 

— 

S 

80 

0.4 

46D 

Characteristic 

Characteristic 

.1574 

— 

S 

80 

0.4 

41D 

Characteristic 

Characteristic 

.1565 

— 

M 

80 

0.4 

46D 

Primitive 

Primitive 

.1589 

.0225 

M 

80 

0.4 

46D 

Characteristic 

Characteristic 

.1558 

— 


Table 1: Shedding frequencies computed with finite-difference code (FD), mixed fi- 

nite-difference-spectral code (M), and fully-spectral code (S). 


These frequencies were based on the computed pressure at a point in the wake 10 di- 
ameters downstream of the cylinder and one diameter above the wake centerline. This 
location corresponds to the principle one used in the experiments of Sreenivasan (1985). It 
should be noted that the frequencies for the finite-difference solution with primitive vari- 
ables at Mach 0.4 differ slightly from those given by Townsend et al. (1987), in which the 
results were incorrect due to a misinterpretation of the output from the FFT subroutine. 

Experimental measurements of the primary shedding frequency in the wake of a circular 
cylinder at low Reynolds numbers were reported by Roshko (1953). An approximation to 
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the data was given by the equation 

Sti = 0.212 fl - J 1 ' 2 -) • ( 3 - 4 ) 

For a Reynolds number of 80, this equation yields fjti^O.1558. As shown in table 1, the 
computed values of St x are very close to the experimental value in all cases. The largest 
deviation is 3.47% and most of the values are within 2% of the experimental frequency. 
Therefore, the primary vortex-shedding frequency in the wake is well-predicted with either 
set of boundary conditions. 

Plots of the time history of the pressure and the resulting frequency spectrum at the 
selected point in the wake flow from the finite-difference calculations are shown in figure 4 
for both sets of boundary conditions. The corresponding plots for the fully-spectral calcu- 
lations with characteristic boundary conditions are shown in figure 5. As shown in figure 4, 
the solution with primitive-variable boundary conditions produces a low frequency which 
is not seen in the spectrum from the calculations with the characteristic boundary condi- 
tions. This frequency is not a subharmonic of the primary shedding frequency and was 
not found in any of the solutions using any of the three codes when characteristic bound- 
ary conditions were employed. Furthermore, the amplitude of the secondary frequency 
remained relatively constant at all locations for which the spectrum was computed, in- 
cluding a location 10 diameters upstream of the cylinder. At this point, the amplitude 
of the primary shedding frequency was significantly lower than in the wake so that the 
secondary frequency was the dominant frequency in the signal. These factors suggest that 
the secondary frequency is a spurious numerical artifact. 

Calculations were also made with the finite-difference code using the primitive- variable 
boundary conditions to determine the effect of Mach number on the computed frequencies. 
Computations were made at Mach numbers of 0.2 and 0.1 for a Reynolds number of 80. 
As shown in table 1, there was only a small effect on St x as Mach number was lowered. 
However, as Mach number was decreased by a factor of 2, Sti increased by a factor of 
approximately 2. 

To further establish that the secondary frequency is a numerical artifact, calculations 
were made for Mach 0.4 flow at a Reynolds number of 20. Under these conditions, the flow 
in the wake should be steady, consisting of two symmetric counterrotating vortices just 
behind the cylinder. Calculations with the finite-difference code produced the expected 


11 


steady wake with both sets of boundary conditions. However, when the primitive-variable 
boundary conditions were used, a secondary frequency appeared in the frequency spectrum 
of the pressure in the wake with a value of 5^=0.0158. This spurious secondary frequency 
did not appear when characteristic boundary conditions were used. 

In summary, the calculations presented in this section demonstrate that a secondary 
frequency can be introduced into the solution by the far-field boundary conditions. By 
coincidence, the value of this frequency was very close to that found experimentally by 
Sreenivasan (1985). This frequency was not found when the boundary conditions are 
properly formulated in terms of characteristic variables. 

4. Analysis 

As shown in §3, the secondary frequency disappeared when boundary conditions based 
on the characteristic variables were used at the outer boundary, demonstrating that the 
secondary frequency computed with the primitive- variable boundary conditions was, in 
fact, of numerical origin. In this section, a theoretical explanation is given for these results 
using a modal analysis of the effect of inflow boundary conditions on the stability of small 
perturbations to the flow. It is shown that the use of boundary conditions based on 
noncharacteristic variables causes a temporally-periodic perturbation. Furthermore, it is 
demonstrated that the use of characteristic variables eliminates this periodic disturbance. 

The theoretical study is based upon a one-dimensional treatment of the Euler equations. 
The use of the inviscid form of the governing equations is justified at distances far from 
the cylinder since the Reynolds numbers under consideration place the flow well outside 
the limits of the Stokes and Oseen regions. The assumption of one-dimensionality is more 
problematical. However, it is a reasonably good representation for the calculations in which 
a finite-difference method was used in the radial direction and the Fourier collocation 
method was used in the circumferential direction. In fact, it will be shown that the 
theoretical predictions of the secondary frequency are in good agreement (within 15%) 
with the computed values for this case. The agreement with the computed results from 
the fully finite-difference code is not as good, although the trend with Mach number is 
well-predicted. 

The analysis considers a region of flow along the x axis sufficiently far upstream from 
the cylinder. The one-dimensional form of equation (2.1), under the assumption of inviscid 
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flow then becomes 


dU 

dt 



(4.1) 


where 


U = 


P 

m , 
E 


F = 


m 

f+P 


[ u(E + p) \ 

and m = pu. In terms of the dependent variables in U, the equation of state, equation 
(2.7), is 

r i „2 l 

(4.2) 


P = (l “ 1) 


1 m 2 
E- -( — ) 

2 k p ’ 


Equation (4.1) can be rewritten in nonconservation form as 

f + ^= 0 ' 


(4.3) 


where A(U) = §£ is the Jacobian of the flux with respect to the solution vector U. 
Linearizing about the steady free-stream conditions, this equation becomes 


l(SV) + A(U a ,)Al6U)=0, 


(4.4) 


where 6U = U - is the perturbation vector. The matrix A^Um) = ^ ias 

three eigenvalues a\ = u <*, — Coo, 02 = tioo + c oo> and a s = u oo • The corresponding 

eigenvectors in terms of the conserved-variable perturbations are 


Ri = -(6m - Uqo 6p) + 

R2 — ( 6 m u 00 6p) ~t~ 


Co o 
7-1 


^ 1 -u^,£p _ «oo 6m + 6E 

, z 

^ulcSp-iiooSm + SE 


and 


R$ — c<x>6p 


7-1 


^ u lo$P - Uoo 6m + 6E \ . 


(4.5) 

(4.6) 

(4.7) 
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Furthermore, using a linearized version of the equation of state, i.e., 

<5p = (7 - 1) \ u lo&P ~ UcoSm + <5.e| , (4.8) 

equations (4.5)-(4.7) may also be written in terms of the primitive-variable perturbations 
(6p,6u,Sp), giving, 


Ri = 

Poo^U “I - 

Coo 

(4.9) 

i? 2 = 

Poo^U “f” ^p, 

Coo 

(4.10) 

Rs = 

Coof>P f>P- 

Coo 

(4.11) 


and 


Any perturbation imposed on the free-stream solution will evolve as a combination of these 
eigenvectors. This evolution can be studied locally near the outer computational boundary 
by employing the modal analysis developed by Gustafsson, Kreiss, and Sundstrom (1972) 
and Osher (1969) to study the stability of numerical approximations to initial-boundary- 
value problems. 

Consider an inflow boundary point with subsonic flow crossing the boundary into the 
solution domain. At this point, a 2 = «oo + Coo > 0 and a 3 = Uoo > 0. Therefore, the 
characteristic variables, iE 2 and i?s, corresponding to these two incoming characteristics 
must be specified. The third one, R x , corresponds to a 1 = — Coo < 0 and thus exits 

the domain. To model the situation near the inflow boundary, consider the governing 
equations in terms of the characteristic variables, viz., 


dR x 


+ («c 


dR x 


Coo) 'dx 


dt 

dR 2 ( x dR 2 

"H (Uqo + Cqo) 


dt 


or 


dR, 




dR 3 
dt 


dR, 


+ u 0 


dx 

dRs 

dx 


0, 

0, 

0, 


= 0, s — 1,2,3. 


(4.12) 

(4.13) 

(4.14) 

(4.15) 


dt ' 3 dx 

In the one-dimensional case, the finite-difference algorithm described in §3 is equivalent 
to a Lax-Wendroff scheme and can be written as 

°‘ A1 (»;?. - «-;?.) + (»;?. - *»;•" + <0 . / > «, (■*.») 


a.n+1 

Wj' = Wj 


2Az 
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where w"' n — w 5 (jAx,nAt) is the finite-difference approximation to R 3 (x,t). Equation 
(4.16) is solved with inflow-type boundary conditions imposed at x b0 undar V — (jAi),=o 
which are of the form 

{R 2 ~cxRi ) j= o = 0, (4.17) 

{Ri-0Ri)i=o = ( 4 - 18 ) 

and 

(.Ri -+- oR 2 + (.R$) j—Q = ( R\ 4* oR 2 + cR^j—i- 

Thus, equation (4.19) is a zeroth-order extrapolation on Ri- 
The modal solution ansatz for this case is 

R„ = A s K,{z n s = 1,2,3. (4.20) 

Substituting this solution into equation (4.16) yields the three modal characteristic equa- 
tions given by 

(A? - AX + 2(1 - z - A*)/c s + (Aj + A,) = 0, s = 1, 2, 3 (4.21) 

where \ 3 = Substituting equation (4.20) into equations (4.17) to (4.19) gives the 

modal representation of the boundary conditions, 

A 2 = aAi, {Art 0) (4.22) 

As = flAi (4.23) 

and 

.Ai(l — /Ci) + (tA 2 {1 — /C 2 ) + cAs{l — /C 3 ) = 0. (4-24) 

Since A x t 0? equation (4.24) may be rewritten as 

(1 — /ci) + acr(l — k 2 ) + /3e{ 1 — /c 3 ) = 0. (4-25) 

It is necessary to find z such that | z |= 1 and /c 4 (s = 1,2,3) such that | /c, |< 1 which 
solve equations (4.21) and (4.25) simultaneously. It should be noted that | z |= 1, z t 1> 
corresponds to a solution that is purely periodic in time. If the search for such solutions fails 


15 


(which it does, as it will be shown, when a = 0 =0, i.e., for characteristic specification), 
then the boundary treatment does not introduce any spurious frequency. If, under a given 
set of boundary conditions, there is a solution to equations (4.21) to (4.25), then the 
phase of z gives the temporal frequency of the numerically-introduced perturbation. Of 
course, it is hoped that at least one of the /c’s will be nearly 1 in magnitude. Otherwise, 
the temporally-periodic perturbation will decay spatially as |/cp. In fact, in all of the 
calculations with the model, it is always the case that .995 <| k .2 |< 1. 

In the actual numerical computations presented in §3, the noncharacteristic boundary 
conditions (translated to the present one-dimensional model) consisted of specifying the 
velocity and temperature ((6u) ;=0 = 0, (<5T) J= 0 = 0) and performing the extrapolation 

Mi-o = (4.26) 

1 oo 

This is equivalent to 

[bp)j=o = 

This is not a “pure” extrapolation on the density. Using (6u) J= o = {6T)j = o = 0 in equations 
(4.17) and (4.18) leads directly to a = 1 and 0 = 7 - 1. If equation (4.27) were a pure 
extrapolation on p, then from the relation 

2coo^p = Ri + R 2 + 2J?3 (4.28) 

(see equations (4.9) to (4.11)) it would follow that a = 1 and e = 2. However, it should 
be noted that while 6u does not appear in equation (4.27), and hence a = 1 in (4.19), 
e must have a value which depends on the solution inside the domain. Therefore, it is 
necessary to seek solutions to the nonlinear system (4.21) and (4.25) for a given A, such 
that | z |= 1,2 ± 1, | k s |< 1, and acr =1. The solution that satisfies z — e' e (0 ^ 0), 

| K a |< l(s = 1,2,3) with ao=l then yields the value of t0 = (7 — l)t. The numerical 
solution of equations (4.21) and (4.25) can be sensitive near | /c, |= 1 so double precision 
was used in the numerical solutions of these equations. For a wide range of M (and hence 
A), it was found that e0 = —1.24 ± 0.2, giving confidence in the model 2 . The A used in 

2 In fact, it can be easily shown that the model boundary conditions (4.17)- (4.19) , with a = a = 1, 
^ = 7 — 1 , and £ found from the solution to (4.21) to (4.25), correspond to a primitive-variable boundary 
condition of the form (p) J= o = {p)j=i^F ± + a(Ar max )‘ 2+m where a = 0(1) and m > 0. Thus, the analysis 
is performed on a boundary condition wliich, within the second-order accuracy of the finite-difference code, 
cannot be numerically distinguished from the one actually imposed i.e., equation (4.26). 


(sr-)(*P)i-i + 

^ OO OO 

(£-)(«*>),=> + £=-(MV.- (4.27) 

00 00 
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the theoretical prediction, A = (u + c)f£, was based on the actual time step used in the 
Navier-Stokes codes and Ax was taken as the radial cell size (Ar maI ) nearest to the inflow 

boundary. 

The present theoretical model most closely corresponds to the Navier-Stokes code with 
the finite-difference method in the radial direction and the Fourier collocation method 
in the circumferential direction. A comparison of the theoretical Strouhal number from 
the one-dimensional model with the one obtained from this Navier-Stokes code is given 
in Table 2. In view of the fact that the analysis is based on a one-dimensional model, 
the agreement between the values of St 2 predicted by the model and the computed values 
is reasonably good. Furthermore, the model predicts the trend of decreasing St 2 with 
increasing domain size. 


N 

max 

^^code 

^ model 

{St^code 

(^ 2 ) model 

100 

.41 

.0053 

.0452 

.0244 

.0296 

100 

.46 

.0053 

.0403 

.0225 

.0264 


Table 2: Comparison of secondary frequencies predicted by the analysis and computed 
using the mixed spectral-finite-difference code. Moo 0.4. 

Finally, it is seen that if the characteristic boundary conditions, (i.e., a = (3 — 0 in 
equations (4.17) to (4.19)) are used, then the only solution to (4.21) through (4.25) is 
A 2 = As = 0, /ci = 1, and z = 1. It can be shown by perturbation analysis that it is not 
a generalized eigenvalue. Thus, any perturbation from the boundary, Ri = AiK.[z n = A x 
will remain constant and small. 

5. Concluding remarks 

A numerical study has been conducted using three different codes to compute the 
unsteady flow about a circular cylinder in low-speed flow at low Reynolds numbers. One 
of these codes used a finite-difference method to solve the two-dimensional time-dependent 
compressible Navier-Stokes equations, the second used spectral techniques, and the third 
used a combination of these two methods. With stable and consistent boundary conditions, 
all of these methods were able to predict accurately the major features of the flow such 
as the vortex-shedding frequency. However, it was found that certain far-field boundary 
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conditions which used extrapolation of the primitive variables introduced an additional 
temporal frequency into the solution. By coincidence, the value of this frequency was 
very close to that found experimentally by Sreenivasan (1985). The use of characteristic 
variables in the far-field boundary conditions eliminated this frequency from the solution. 
An explanation of this behavior was provided using an analysis based on a simplified model 
for the boundary conditions. These results illustrate the fact that great care must be 
taken in interpreting the results from numerical simulations, and that while the secondary 
frequencies are spurious in the infinite-domain case, the possibility exists that they may 
be found in a confined space such as a wind tunnel. In fact, Sreenivasan 3 found that 
the secondary frequency disappeared when the tests were performed in a wind tunnel 
where the distance from the cylinder to the upper and lower tunnel walls was 250 cylinder 
diameters. This distance is about seven times that used in the original tests. Insight into 
the presence of the multiple frequencies in the original experimental results would require 
three-dimensional computations which included the wind-tunnel walls. 
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(a.) Primitive-variable boundary conditions. 



(b.) Characteristic-variable boundary conditions. 


FIGURE 2. Lines of constant pressure (-^). Finite-diflerence code. Contour increment 
=0.01. M x = 0.1, Re „ iD = 80. 
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FIGURE 3. Lines of constant pressure (-£-). FuIIy-spectral code. Contour increment 
=0.01. Moo = 0.4, Re^'O = 80. 
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step number 

Primitive-variable boundary conditions. 



step number 

Characteristic-variable boundary conditions. 


(a.) Time history of pressure. 
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(b.) Frequency spectrum of pressure. 


FIGURE 4. Effect of boundary conditions on pressure in wake region at location 10 
cylinder diameters downstream of cylinder and one diameter above wake centerline. Finite- 
difference code, psd = power spectral density. 
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step number 

(a.) Time history of pressure. 



(b.) Frequency spectrum of pressure. 


FIGURE 5. Computed pressure in wake region at location 10 cylinder diameters down- 
stream of cylinder and one diameter above wake centerline. Fully-spectral code, psd = power 
spectral density. 
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